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Abstract 

An asymptotic analysis of the Gunn effect in two-dimensional samples of bulk n-GaAs with 
circular contacts is presented. A moving pulse far from contacts is approximated by a moving free 
boundary separating regions where the electric potential solves a Laplace equation with subsidiary 
boundary conditions. The dynamical condition for the motion of the free boundary is a Hamilton- 
Jacobi equation. We obtain the exact solution of the free boundary problem (FBP) in simple 
one-dimensional and axisymmetric geometries. The solution of the FBP is obtained numerically 
in the general case and compared with the numerical solution of the full system of equations. The 
agreement is excellent so that the FBP can be adopted as the basis for an asymptotic study of the 
multi-dimensional Gunn effect. 

PACS numbers: 73.50.Fq, 73.61.Ey 
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I. INTRODUCTION 



Excitable media exhibit a large response to a sufficiently strong disturbance from their 
only stable stationary homogeneous state. This feature makes them ideally suited to sustain 
propagation of pulses or wave trains Ij . Examples are the propagation of an action potential 



along the axon of a nerve , the propagation of a grass fire on a prairie, pulse propagation 
through cardiac cells react ion- diffusion or ecological systems Semiconductor 
systems displaying negative differential resistivity in their current-field characteristics are 
also excitable systems albeit they have peculiar features due to the long range character of 
the electromagnetic interaction . Thus dc voltage bias conditions lead to pulse recycling 
(at contacts) and motion that give rise to self-sustained oscillations of the electric current, 
the so-called Gunn effect, jsl- While most of the theoretical and experimental studies of these 
phenomena deal with one dimensional geometries of samples with attached planar contacts, 
recent experiments and numerical studies 0, 0] have considered rectangular samples with 
point contacts. In this case, many unusual oscillatory patterns are found [l,!^- 

A large part of the literature on pulse propagation is devoted to the mathematical descrip- 
tion of their motion on one-dimensional unbounded domains. In the case of self-oscillations 
in semiconductor systems, such a description is the basis of asymptotic analyses of pulse re- 
cycling and motion, both in one-dimensional (ID) 10| or axisymmetric two-dimensional 
(2D) samples 0, llll- These studies exploit that the electric field has only one relevant 
component whose integral yields the voltage difference between contacts. The situation is 
very different in more general geometries and new ideas need to be brought in. In this 
paper, we reduce pulse propagation (far from contacts) to the motion of a free boundary 
(FB) separating regions where the electric potential is a harmonic function. The FB obeys 
a Hamilton- Jacobi equation (HJE). On the FB, continuity and jump conditions hold, and 
additional conditions on contacts and sample boundaries are needed for the problem of find- 
ing the FB (free boundary problem or FBP) to have a unique solution. On simple ID and 
axisymmetric geometries, the HJE can be solved exacly. Its solution describes very well the 
motion of a discontinuity of the electric potential representing the pulse far from the bound- 
aries, as comparison with the numerical solution of the full system of differential equations 
shows. This is also true of the general 2D case, but now the solution of the FBP has to be 
obtained numerically. In all cases, recycling and annihilation of pulses at contacts have to 
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be described separately from the FB motion. 

The rest of the paper is organized as follows. In Section m we present the governing 
equations of the Kroemer model for the Gunn effect in two-dimensional samples of n-GaAs. 
An asymptotic derivation of the FBP is given in Section 11111 Section IIVI contains the exact 
solutions of the FBP in the ID and axisymmetric cases. Numerical solutions of the FBP 
in the general 2D case and comparisons with the numerical solution of the full system of 
equations are presented in Section |Vl The last Section contains our conclusions. 



II. EQUATIONS AND BOUNDARY CONDITIONS 



The Kroemer model 12| consists of the following equations and boundary conditions (in 
dimensionless units) for the concentration of free carriers (electrons), ra, and the electric 
potential, y?: 

dn 



dt 



+ V ■ {nv-5Vn) 



n 



v{E) = E 



1 + VsE^ 



l + E^ 
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Here (P) and Q are the charge continuity and Poisson equations, respectively. The dimen- 
sionless electric field is E = Vy? and E = \E\. In these equations, the electron density 
has been scaled with the uniform concentration of donor impurities in the semiconductor, 
Nd = 10^^ cm~^, and the electric field with the field characterizing the intervalley transfer 
responsible for the negative differential mobility involved in the Gunn oscillation, E^i = 3.1 
kV/cm. Distances and times have been measured with the dielectric length and the dielec- 
tric relaxation time, h = eEji/ {eNo) ^ 0.276/im, Ii/{pqEr) ^ 1.02 ps, respectively (/io is 
the zero-field electron mobility; see, e.g., lOl] for details). The unit of electric potential is 
Ejili ^ 0.011 V. The carrier drift velocity of Eq. (jH}, v{E), is already written in dimen- 
sionless units, and it has been depicted in Fig. 1 of Ref. We assume that the diffusion 
coefficient is constant, 6 ~ 0.013 (at 20K). In the rest of the paper we assume also a zero 
saturation velocity: Vg = 0. 



Boundary and bias conditions need to be imposed at the interfaces between semiconduc- 
tor and contacts, Sc,a, and on the outer boundary of the semiconductor boundary Sq. Our 
boundary conditions (j3]) and © assume that the normal components of electron current 
density and electric field are proportional at the semiconductor-contact boundary (Ohm's 
law) y, (in these equations, N is the unit normal to Sc,a; directed towards the semicon- 
ductor). For simplicity, we choose all contact resistivities p to be equal. Bias conditions are 
chosen to be = at the cathode Sc (injecting contact) and ^9 = $ (the applied voltage) 
at the anode (receiving contact). If part of the semiconductor boundary does not have 
attached contacts, the corresponding boundary conditions are zero flux ones, as in Eq. 0. 
Typically 5 > is very small, so that diffusion matters only inside boundary layers near the 
contacts or inside thin shock waves 0, 2. The latter are charge accumulations that will be 
treated simply as discontinuities of the electric field j^. Thus diffusion effects may be left 
out of the conservation equation when interpreting the results. If we set 5 = 0, the first 
boundary condition in Eq. and the second one in Eq. ^ should be omitted. 

We can write an Ampere's equation for the total current density (electronic plus displace- 
ment), j, by eliminating n from using 

V ■ J = 0, with 
^ ^ BE 

j = (l + VV)t'-5V(VV) + ^. (7) 

III. DERIVATION OF THE FREE BOUNDARY PROBLEM 

Let us consider a rectangular sample with circular contacts whose radii Vc are large but 
much smaller than the distance between contacts, 1 < -C L. The current density varies 
slowly and follows adiabatically the electric field profiles in the semiconductor except during 
brief periods in which new pulses are shed from the cathodes. Close to a cathode located 
at the origin, the electric field and the current density are approximately axisymmetric and 
we can use the results of Ref. [ill. j = Jr/r^, r = \r\, E = Ei{J/r)r/r. Ei{j) and E2{j), 
with El < E2, are the two positive zeros of the function v{E) — j, with v{E) = \v{E)\. The 
maximum value of |j| during self-oscillations is somewhat larger than = 0(1) at which 
E2{j) = pj- Correspondingly, the maximum value of J is Jc = jc^^c = 0{rc) and far from the 
cathode, r ^ r^, J <^ r holds. This means that E ~ Ei{J/r) ^ J/r <^ 1 and v{E) ^ E. 
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When Vg = the pulses move slowly over large regions of the sample in which the field is 
stationary and small: E <^ 1. Notice that v{E) = E — E^/{1 + E^), which implies v{E) 
to be approximately linear on a wide range of field values, E^ <^ 1. We conclude that 
v{E) ^ E except near the contacts and inside pulses. In these outer regions, space and time 
derivatives can be neglected in Eq. ((7j), which implies j ~ v{E) ^ E there. Thus divj = 
yields V'^f = and the electric potential ip is a harmonic function outside pulses and contact 
regions: 

VV = 0. (8) 

Let us now consider the pulse interior. A pulse is a narrow region of high electric field 
bounded by a leading front and a trailing front which is a shock wave. Outside the pulse 

^ 1 as explained before. The leading front is a region at which n = l + V- -E'~0. Since 
we are describing the pulse far from the contacts, r ^ Tc ^ 1, the electric field is essentially 
normal to the pulse, A^. Then Ej^j = E ■ N ^ ryj{t) — r , where r measures displacement 
along the normal to the front and rw{t) yields the front location. The velocity of the leading 
front is dr^/dt = = j ■ N, according to Eq. (jZj). The back of a triangular pulse of height 
(E^ — E^) (the trailing front) is a shock wave with speed given by the equal area rule [uj 

V{E,, E.) = ^-1^ v{E) dE ~ ^, (9) 

where we have used that E^ ~ Ei 13] as ^ 1. Then the trailing front velocity is small 
and small waves move faster than large ones. A key observation is that the pulse is narrow 
and it can be substituted by a curve on a length scale of the order of the distance between 
contacts, L. This is clear if leading and trailing fronts of the pulse are circular ^|. Then 
the bias $ = 0{L) is the integral of the electric field from the cathode to the anode and 
the pulse width (equal to its height) is {E^ — E_) = 0(-\/$) ^ In the general case, the 
pulses are circular during a large part of their lives [if and we shall assume that their widths 
remain much smaller than L even when their shapes are no longer circular. Then we assume 
that the pulses are curves V given by the equation: 

Vr(f,t) = 0. (10) 

Clearly, there is a finite voltage drop across the pulse, ~ E\l2 = 0($), which means 
that the electric potential has a jump discontinuity at F: 

^ = (11) 



Here (/?_ and (p+ are the limiting values of as x approaches F from the region inside or 
outside r, respectively. The relations divj = and j ~ -E imply that the normal component 
of the electric field (and therefore the normal derivative of the electric potential) is continuous 
across F: 

j^ = {N-Vip)^ = {N-Vip).. (12) 

This jN is also the velocity of the leading front of the pulse along its normal, which 
is nearly equal to that of the trailing front, V given by Eq. during most of the pulse 
lifetime. The pulse velocity can also be obtained by differentiating Eq. (|Tnjl with respect to 
time: 

dW -L^^^ dx 

Since N = VW/\VW\, the normal component of the pulse velocity, Jn, is 

(13) 

dt \VW\ dt 

Using Eqs. Jn = V and (fT^. we obtain the following equation for the position of the 
FB F: 

' ' on W{x,t) = 0. (14) 



Thus we have posed the following FBP: 

The electric potential (f{x,t) is a harmonic function inside and outside the FB T, with 
boundary conditions (0) and ^ on the semiconductor boundaries. On the FB T, im- 
plicitly given by W{x,t) = 0, (f has a jump discontinuity [(p] and its normal derivative 
satisfies 

(iV-Vvp)+ = ^= = (iV-V^)_, 



where N = 'VW/\'VW\. Furthermore, the FB obeys the following HJE: 

dW vrlWl 



on W = 0. 



The conditions on the normal derivative of the electric potential at the FB are equivalent 



to: 



7l\VW\ 

4,/2M 



onW = Q. 

The HJE ()14|) can be solved by the method of characteristics (the Hamilton equations). 
To derive them, we just take a partial derivative of the HJE with respect to and a partial 
derivative with respect to y. The results are 

d dW TT fdWd'^W dW d'^W\ 

dt dx A\VW\sj2[ip\ dx"^ 9y dxdy ) 

and a similar equation for dW/dy. The corresponding characteristic equations for these 
first-order quasilinear partial differential equations for p = dW/ dx and q = dW/ dy are 



dx 4a/2[^ 



dt y/p'^ + q'^ 



dy _ 4^/2M 
dt ~ Vl^T^"^ 

dp TT d[ip] 

dq TT d[{p] 

0. 



dt 

In these equations s is arc length on the FB F, and we have used that d[(p]/dx = 
— g((9[(/9]/(9s)/A/p^~+~g^ and d[(p]/dy = p{d[{p]/ds)/ y/p'^ + q^ on F. These expressions can 
be straightforwardly derived by using a local coordinate system on F with basis vectors 
N = VW/\VW\ and f = {-dW/dy,dW/dx)/\VW\. The jump [if] depends only on the 
arc length on F and t because it is defined only for (x, y) E T [these (x, y) G F have zero 
projection onto A^]. The last equation for W follows from the chain rule, the Hamilton 
equations for x and y and the HJE: 

dW dW dWdx dW dy 
= \ \ ^ 

dt dt dx dt dy dt 
dW ttIVWI ^ 
dt 

The characteristic equations can be used to find W{x, y, t) given an initial condition 
W{xo,yo,0) = Wo{xo,yo) such that the FB is described initially by Wo{xo,yo) = 0. Let 
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us assume that [ip] is a known function of s and t. In principle, we can find the solutions 
of the above equations with initial data x = xq, y = yo, p = dWo/dxo and q = dWo/dyo- 
The result is a two-parameter family of solutions x = X(t;xQ,yQ), y = Y(t;xo,yo). Let us 
assume that we can invert this transformation for each t > (which should be true for t 
sufficiently small), Xq = C,{x,y,t), y^ = ri{x,y,t). The solution of the HJE is W{x,y,t) = 
Wo{^{x,y,t),ri{x,y,t)) because W is constant over the characteristics. Once W is found 
for a given the Laplace equation can be solved for the electric potential in the different 
regions of the sample separated by the FB W = 0. Inserting these solutions in the definition 
of the jump [ip], we find an equation for this jump. It seems clear that we can implement 
this procedure numerically to device an explicit method such that W and [(f] are calculated 
at time t + At knowing their values at time t. In particular, we do not need to find (xq, yo) 
in terms of {x,y,t). All we need is to know the instantaneous location of the FB W = 0, 
thus we only need to know the evolution of those {xo,yo) that are on Wq = 0. For each 
t > 0, the locus of such x = X{t;xo,yo), y = Y{t;xo,yo) constitutes the FB. More details 
on the numerical implementation of these ideas are given in Section |^ 

The FBP describes the motion of a pulse far from contacts and other boundaries or 
pulses. To obtain a complete asymptotic description of the Gunn self-oscillations, we have to 
supplement its solution with a local description of the field near the contacts and boundaries 
and a description of pulse collisions. In particular, new pulses are shed from the cathodes 
as the normal component of the current densitv there surpasses a critical value jc which is 

. There are cases in which two pulses collide and 
In these cases our construction of the moving FB 
breaks down. What do we do then? Consider for instance two circular pulses that become 
tangent at a point {xi,yi) at time ti > 0. Clearly there are two different initial points 
(xo,i/o) that have evolved towards {xi,yi) and W{x,y,t) is no longer univalued. Numerical 
simulations of the complete system of equations show that the two pulses merge and adopt 
an eight-shaped form; see Fig. 7 of Ref. ui To mimic this situation with our FBP, we should 
stop the simulations and start with a new FB at t = ti that is an eight-shaped simple curve 
with a hole at the tangent point of the two old pulses. The new FBP should now have a 
unique solution. 



the same as in the axisymmetric case [1 
merge and cases in which a pulse splits jj 
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IV. EXACT SOLUTIONS OF THE FREE BOUNDARY PROBLEM IN SIMPLE 
GEOMETRIES 



There are two simple geometries in which the FBP can be solved exactly: parallel planar 
contacts attached at the ends of a rectangular sample (ID case) and the Corbino geometry 
of two concentric circular contacts with the sample in between (axisymmetric case). Let us 
call region A that comprising the cathode and region B that comprising the anode. 



A. ID geometry 

Then the electric potential depends only on the coordinate x, the cathode is located at 
X = 0, the anode at x = L and the FB is a moving point Xs{t). The electric potential obeys 

m (0,Xs), V5yi(0,t) = 0, — — (a;s,t) 



ox-^ ox 



m [Xs, L), lpb[L, t) = $, -7r—{Xs, t) 



dx^ ' ' ' ' dx ' - ' 4y'2[v?] 

The solutions are 

4./2M 



(Pb{x, t) = — ^= (x - L) + $. 



The jump in the potential, [(p] = ipB{xs,t) — ipA_{xs,t) is independent of t and Xs, and it 
solves the following equation: 



$ ; L. 



4V2M 

Setting a = \f\(p\ and = we obtain 

«'=(0«-^)^- (15) 

Depending on the values of and L this equation may have zero, one or two positive 
solutions. If there are two solutions, an argument due to Volkov and Kogan 14] shows that 
the pulse with smaller [ip] is unstable. The FB Xs{t) can be found by solving the dynamical 
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FIG. 1: The solid lines indicate the electric potential and field of an advancing ID pulse (far from 
the contacts) calculated by numerically solving the Kroemer model. They agree very well with 
the approximations ipA{x,t) and ipB{x,t) (dashed lines). We have used the nondimensional units 
defined in the text. 



HJE: 



dW 
'~dt 



IT 



dW 

dx 



V2b] 

Let us assume that the initial profile W{x,0) = Wo{x) is monotone increasing and that 
it vanishes at a position Xs{0) G (0,L) corresponding to the pulse location at time t = 0. 
For small enough t, we then have dW/dx > and we can ignore the absolute value in the 
previous equation. Its solution is then 



4./2 



W{x,t) = Wo ^ 

Notice that we have dW/dx > for all t > 0. Since W{xs,t) = 0, the previous solution 
yields 

vr 



Xs{t) = Xs{0) + 



4./2 



(16) 

Fig. [T] compares (pA{x,t) and ipB{x,t) to the electric potential of an advancing pulse 
calculated by numerically solving the exact system of equations. 

The FBP has yielded the same approximation to the complete ID problem as indicated 
in Ref. for the motion of a pulse far from the boundaries. When the pulse arrives at the 
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anode x = L, it starts disappearing there and the current density increases until it surpasses 
jc- Then a new pulse is shed at x = 0; see Ref. for details. 

B. Corbino geometry (axisymmetric case) 

The potential depends only on the radius r measured from the center of the cathode. 
Solving the Laplace equation d[r dip/dr]/dr = at both sides of the moving pulse of radius 
rs{t), we find 



(pA{r,t) = ^Zl- log (-) 



Mr,t) = log ) + $. 



c 



The jump in the electric potential at r<j is now given by the following equation: 

M = '^ r== log 

or equivalently 

«3 = _ ^ log ("Jill ) r,. (17) 



4^2 



for a = y . Notice that explicitly appears in these equations and that [yj] decreases as 
the pulse advances (and therefore increases); cf. Ref. [ill The HJE can be solved as in 
the ID case and its solution yields 



r,(t) = r,(0) + ^/^[^]'^rft. (18) 

In this case, Eqs. (fTTj) and (fTHj) for [ip] and Tgit) need to be solved simultaneously. 

The stage of a self-oscillation described by the previous FBP corresponds to having a 
single pulse far from the contacts. See Ref. |ll| for a fuller description of self-oscillations in 
this case. 

V. NUMERICAL RESULTS 

To test our FBP formulation, we shall consider the relatively complicated geometry of 
Fig. 7 in Ref. Q (reproduced here as Fig. |21 to facilitate comparison with the results of 
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FIG. 2: Density plots of the solution of the Kroemer's model (with Vs = 0) in a square of side 
1=20 with four circular contacts forming the vertices of a square of side d = 4 located at the center 
of the sample. Cathodes have potential 99 = and anodes have <p = 10. Our dimensionless units 
have been defined in Section Ull 



numerically solving the FBP) corresponding to Vs = 0. The sample is a square of side / = 20 
with two cathodes at potential (f = and two anodes with (f = 10. The circular contacts 
(of radii 0.5) are at the vertices of a square of side d = 4 located at the center of the sample. 
Then the separation between contacts is L = 3 and the distance from contacts to the border 
of the sample is 7.5. Notice that dipole waves are emitted from the cathodes. Immediately 
after their emission, the waves are circular. As they approach each other, the waves become 
elongated and merge forming an eight-shaped connected curve that grows until it reaches 
the anodes. 



A. Free boundary problem 

We shall now explain the results obtained by solving numerically the FBP. Details of 
the numerical method will be given later. Fig. 01 shows the evolution of the FB separating 
the two regions of the sample, inside and outside the boundary. Notice that the numerical 
solution of the FBP closely resembles the numerical solution of the full Kroemer model 
depicted in Fig. |21 In the two first frames of Fig. ^ the FB consists of two circumferences 
corresponding to the dipole waves nucleated at the cathodes. In the third frame, the curves 
collide and then merge forming an eight-shaped closed curve as shown in the remaining 
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FIG. 3: Time evolution of the FB (black curve) separating the two regions of the sample, inside 
(clear grey) and outside (dark grey) the boundary. The anodes appear in white. 
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FIG. 4: 3D plot of the electric potential surfaces (pA{x,y,t) (lower surface, inside the FB) and 
fB{x,y,t) (upper surface, outside the FB) at the time corresponding to the last frame of Fig. El 
Our dimensionless units have been defined in Section Hll 

frames of Fig. |21 Fig. H] shows the electric potential distribution in each region (inside and 
outside the FB) corresponding to the last frame of Fig. El 
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FIG. 5: Dimensionless velocity of each point of the FB at the dimensionless time corresponding to 
the last frame in Fig. |31 and calculated from the electric potential distribution showed in Fig. ^ 

By using Eq. we see that each point of the FB moves with velocity 



Fig. El depicts the velocity of the points at the FB in the last frame of Fig. El The curve 
is symmetric and Fig. El shows that the FB moves faster at the points located in the left- 
upper and right-lower corners of the sample, in agreement with the numerical solution of 
the Kroemer model. 

Let ti be the time at which two dipole waves created at the cathodes touch at a point 
(as in the third frame of Fig. E}, counted from the time at which dipole waves are emitted 
at the cathodes {t = 0). The velocity of the points at the FB is shown at three different 
times in Figures El (0 < t < ti) and [71 {t > ti). Notice that the velocity of the points near 
the center of the sample in Fig. El is larger than in neighboring points, which explains the 
elongated form of the dipole waves in the numerical solution of the Kroemer model (see the 
third and fourth images of Fig. E}. In Fig. Owe observe that the largest velocity is reached 
at the outer points of the single FB, also in agreement with the numerical solution of the 
full model equations. 
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FIG. 6: Time evolution (from bottom to top) of the velocity of the FB P when t < ti, where the 
topology is composed by three domains and F is made of two circumferences. Our dimensionless 
units have been defined in Section ITll 
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FIG. 7: Time evolution (from bottom to top) of the velocity of the FB F when t > ti. The topology 
is now composed by two domains. 

B. Numerical solution of the free boundary problem 

To solve numerically the FBP, we should solve the PDE governing the time evolution 
of the FB, taking into account that the velocity thereof is determined by the solution of 
Laplace's equation with Neumann boundary conditions on the FB and Dirichlet boundary 
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conditions at the contacts (the electric potential problem, or, briefly, the EPP). 

At each time step, the FB advances at a constant velocity for a short distance from its 
previous position. (Thus we ignore the velocity variation during the short time interval 
between ti and tj + At). At time tj + At, we solve the EPP in the different domains resulting 
from the new location of the FB. This yields the electric potential distribution that is used 
to calculate the velocity of the FB at the next time step. 

The time evolution of the free boundary is calculated by using the so-called fast marching 



method (a special case of the method of level sets). This met 



lod was introduced by Sethian 



in 1996 [1^ and used in a wide variety of applications [l^, Level sets methods 

are very efficient for solving complex problems of evolving interfaces whose topology may 
change. If the velocity of the interface does not change sign, the fast marching method is a 
very fast algorithm indeed. 

The general version of the method of level sets consists of solving the evolution equation 

dW 

-— + F\VW\=0, (19) 

where W{x,t) is a function such that Vr=0 describes the free boundary moving at velocity 
F] cf. Eq. ()14|) . When the sign of F does not change, the FB either expands or contracts 
uniformly as time elapses. In our case, the FB moves away from the cathodes. Then the zero- 
level set iy=0 comprises the points farthest from the cathodes that have been traversed once 
by the FB at a given instant of time. Then we can define an arrival time function T in the 
whole sample: T{x) is the time it takes the FB to arrive at the point x starting from a given 
initial configuration. To find an equation for T, we take the gradient of W{x,T{x)) = 0, 
WW + WtWT = 0, and use Eq. (0 to obtain 

VW - F\VW\VT = 0. (20) 

— > — » 

This equation implies that and VT are colinear vectors and their lengths are related 
by |VVr| = -F|Viy| |VT|. Then we obtain the following Eikonal equation for T, 



1 4. /2 [if] 

-m^^- '''' 

The velocity F as a function of x is evaluated at time t. Once the solution of Eq. (|2H) is 
known at a narrow band about the instantaneous location of the FB at time t, the location 
thereof at time t + At is found by solving T{x) = t + At. 
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FIG. 8: The FB comprises two separate curves defining three regions. 



The fast marching method consists of solving numerically this equation by using upwind 



This choice ensures that the information is always taken from where the solution is already 
known. The fast marching method is consistent with the Huygens principle even when two 
waves collide and adopt an eight-shaped curve as in Fig. |21 or with even more complex 
topologies. The EPP is solved by using an integral equation method based upon Green's 
formula. This yields the solution within a region for a given value of its normal derivative 
at each point of the boundary. To make sure that the nonlinear boundary conditions at the 
FB hold, we implement an iterative process. 

We shall start our simulation from an initial configuration as depicted in Fig. |H1 There 
two waves have been nucleated at the cathodes and have reached their typical circular form. 
The FB consists of two circumferences that divide the sample in three regions, Ai, A2 and B, 
in which we should simultaneously solve the EPP. Implementing the fast marching method, 
we see two waves growing from the initial circumferences until a time ti, when they meet 
at the center of the sample. Then the FB is a connected curve and we have the situation 
depicted in Fig. IHl where there are only two regions A and B. The algorithm detects the 
time ti, adapts itself immediately to the new configuration similar to Fig. Inland it continues 
solving the FBP. 

The accuracy and convergence of the method have been successfully checked by decreas- 



finite differences to approximate |VT|. In particular, we have used the Godunov scheme 




(22) 
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FIG. 9: The FB is a single curve defining two regions. 



ing the mesh size. The computational cost of the method is very low compared to the 
computational and memory effort required by the resolution of the full Kroemer model. 
The order-one fast marching method solves the Eikonal equation in the whole sample with 
0{N log N) operations, where N is the size of the mesh, but we only need to solve the 
Eikonal equation in a narrow band ahead of the FB at each time step. On the other hand, 
the EPF solver carries out 0{N^ + M) operations, where M is the number of points defining 
the FB (at most of order N). 



VI. CONCLUSIONS 



We have studied Gunn oscillations in 2D rectangular samples of n-GaAs with circular 
contacts by solving the Kroemer drift-diffusion model with appropriate boundary and initial 
conditions. By using singular perturbation methods, the motion of dipole waves in semicon- 
ductor samples has been reduced to solving a free boundary problem. Exact solutions of this 
problem have been found in simple ID and axisymmetrical (Corbino) geometries. In the 
general case, the free boundary is numerically found by means of the fast marching method 
which is a special case of the method of level sets. The great reduction in computational 
cost allowed by using this method as an alternative to solving the full Kroemer model would 
enable us to study much larger samples. 
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